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Abstract 

In the present work, we generalize earlier considerations for intrinsic localized modes 
consisting of a few excited sites, as developed in the one-component discrete nonlin- 
ear Schrodinger equation model, to the case of two-component systems. We consider 
all the different combinations of "up" (zero phase) and "down" (vr phase) site ex- 
citations and are able to compute not only the corresponding existence curves, but 
also the eigenvalue dependences of the small eigenvalues potentially responsible for 
instabilities, as a function of the nonlinear parameters of the model representing the 
self/cross phase modulation in optics and the scattering length ratios in the case 
of matter waves in optical lattices. We corroborate these analytical predictions by 
means of direct numerical computations. We infer that all the modes which bear 
two adjacent nodes with the same phase are unstable in the two component case 
and the only solutions that may be linear stable are ones where each set of adjacent 
nodes, in each component is out of phase. 



1 Introduction 



The theme of nonlinear dynamical lattices and especially of the prototypical 
equation of the Discrete Nonlinear Schrodinger (DNLS) type [15] has been 
a widely studied one over the last decade, perhaps especially so due to its 
applications in both the field of nonlinear optics, as well as in that of atomic 
physics. 

The realization of this type of models in the context of matter waves arose 
due to the examination of the dynamics of Bose- Einstein condensates (BECs) 
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in periodic, so-called optical lattice, potentials. In that realm, the reduction 
of the prototypical mean-field model of the Gross-Pitaevskii equation with a 
periodic potential, in the superfiuid regime, via a tight binding approxima- 
tion, naturally yields the DNLS as the canonical lattice model for the BEG 
"droplets" [3]fellT6] . 

On the other hand, the nonlinear optics realization in fabricated AlGaAs 
waveguide arrays [8] has been one that preceded the atomic physics one and 
yielded numerous insights about the interplay of discrete lattice dynamics with 
the effects of nonlinearity, including but not limited to Peierls-Nabarro barriers 
and discrete solitons, modulational instabilities and the lattice modifies them, 
gap solitons, diffraction and diffraction management; see e.g. [2pO|9| ?JT8] for a 
number of relevant reviews. Another area of optical applications which has seen 
considerable development over the past decade is that of optically induced lat- 
tices in photorefractive crystals such as Strontium-Barium- Niobate [7|13fl2] . 
Although the latter setting is not directly modeled by the DNLS (chiefly due to 
the photorefractive, saturable nature of the nonlinearity), it has enabled the 
experimental observation of numerous solitary/localized patterns that were 
predicted in the context of the DNLS; see e.g. the reviews [T5][TT| . 

A topic that has been of particular interest in both the atomic and the op- 
tical applications has been that of multi-component dynamics. Among the 
many recent experimental results related to the interaction of two different 
frequencies or polarizations in optics or of two hyperflne atomic states or two 
different atoms in BEG, we highlight only a few. The flrst experimental obser- 
vations of discrete vector solitons in waveguide arrays were reported in [20] . 
while the emergence of multipole patterns in two-component settings within 
photorefractive crystals was presented in P], and more recently more complex 
pattenrs including even dark-bright lattice solitons were reported experimen- 
tally [B] . In the BEG case, experiments with mixtures of different spin states of 
^^Rb [2BPT] or ^^Na [2H] and even ones of different atomic species such as "^^K- 
^^Rb [22] and ^Li-^^^Gs [25] have been realized and although localized modes 
in optical lattices have not been probed in the form considered herein, certainly 
such experiments appear to be within the current experimental reach [23] . 

Our scope in the present communication is to explore systematically the ex- 
istence and stability properties of two-component localized modes, motivated 
by the above applications. The most fundamental excitations (single site) in 
the coupled DNLS equations model were explored a considerable while ago not 
only in one- but also in higher dimensions; see for a relevant discussion [H] 
and references therein. In fact, in higher dimensions also vortex states were 
considered [17]. However, in the one-dimensional realm, we are not aware of a 
mult i- component generalization of the theory of [27] (for the existence theory, 
see also [1]) which enables the systematic quantiflcation of the existence and 
stability properties of localized modes consisting of an arbitrary number of 
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sites. Here, we present the theory of how to address the problem of existence 
and stabihty in the two-component setting, which paves the way also towards 
the generalization to a higher number of components. Aside from providing the 
general theory, we consider numerous special case examples involving excita- 
tion of two- and also of three-sites in both components with either zero phase 
(an "up" site) or with vr phase (a "down" site). For all of these configurations 
we examine their stability properties as a function, e.g., of their nonlinear co- 
efficients (dubbed Qij), which stand for the self- and cross-phase modulation 
parameters in optics and for the scattering lengths (or more accurately their 
ratios) in the realm of BEC. Interestingly, we find that independently of the 
nature of the second component and of the inter- and intra-component cou- 
pling strengths, should two adjacent sites in any one component be of the same 
phase, the configuration is found to be unstable. The only potentially stable 
configurations consist of nodes with alternating adjacent phases, similarly to 
what is know for the one-component problem [27] . 

Our presentation is structured as follows. Section II presents the general ex- 
istence and stability theory. Section III examines special case examples of 
configurations. Finally, section IV summarizes our findings and presents our 
conclusions, as well as topics for potential further studies. 



2 General theory 



Our fundamental model will be that of the two coupled DNLS equations (see 
e.g. [11127115] ) of the form: 



« Un +{9nW + 9l2\Vn\^)Un + eAUn = 

i Vn +{9l2\Un? + g22W)Vn + eAK = 0. 

The localized modes (often referred to also as "discrete solitons") are given 
by the time-periodic solutions of the DNLS equation 



Unit) = 0„e^('^-2e)t+^Co 

where /i G R, 0n,V^n £ C, n G Z, ^o,Co ^ R is a parameter and {^,(j)n,4'n) 
solve the nonlinear difference equations on n G Z : 
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(/i - gi2\4>n\^ - g22\i^n\^)i^n = e(V^n+l + i'n-l)- 



(3) 



Notice that here we start by examining the simpler case where the sohton pa- 
rameter n (representing the chemical potential in the atomic physics case and 
the propagation constant in the optical case) is identical in the two compo- 
nents. Nevertheless, this is not necessarily the case given the nonlinear nature 
of the coupling and the presence of an invariance under phase transformations 
in each of the involved components. We will comment on the case with distinct 
values of /i at the final section of this presentation. 

Our principal existence result can be formulated as follows. 

There exists eo > such that the solution 0„ and ipn, n E Z, are represented 
by the convergent power series for < e < eo: 



fc=i 



and 



where 



k=l 



(4) 
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(6) 



I (1 - 9lMn\^ - ^12|V^nn0n = 
[ (1 - fi'l2|0nP - fi'22|^nnV^n = 0. 

with ^n, Cn £ {0, tt} for n & S = = = {1,2, ■ ■ ■ , N} which is a finite set 
of nodes of the lattice. 

The above analytical expressions fully describe the anti-continuum limit of 
uncoupled sites (i.e., e = 0). In particular, they correspond to the cases of 
both fields being excited at each lattice node, or of those of one field being 
excited only thereon. The latter case was considered previously in [5], hence 
in the present work, we will focus on genuine two-component excitations. No- 
tice also that in the above considerations (and in what follows), we have set. 
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without loss of generality, /i = 1. Lastly, we should point out that for the 
two-excited-component solutions to be meaningful, it should be the case that 
{912 - fl'22)(fl'?2 - fi'iifl'22) > and similarly that {gi2 - fi'ii)(5'?2 - fl'115'22) > 0. 

The spectral stability of discrete solitons in studied with the standard lin- 
earization 

where \,an,bn,Cn,dn solve the linear eigenvalue problem on n G Z: 



anCin - gii4>lbn " fi'i20n^n(c„ + dn) - e(a„+i + a„_i) = iAa„ 

OinK - gil4>n(^n " 5'l20nV'n(c„ + dn) - e(6„+i + = -iXh^ 

/3nCn - 5'22^„C„ - 5'l20n^n(an + " e(c„+i + C„_i) = iAc„ 
/3ndn - g22i'ldn - fi'l20nV'n(a„ + &„) - e(rf„+i + dn-l) = -iXdn, 

where a„ = 1 - 2gn(l)l - gu^l, /3„ = 1 - 2^(22^/'^ - fi'i20^- 

Let = /^(z, C) be the Hilbert space of square-summable bi- infinite complex- 
valued sequences {un}nez , equipped with the inner product and norm 

(U, V)n = J2 "^nVn, II U ||^= ^ |m„P < OO. (10) 

n6Z nez 

We use bold notation u for an infinite-dimensional vector in VL that consists 
of components Un for all G Z. 

The stability problem dH]) is transformed through the substitution 

an = Un + iWn, bn = Un - iWn, Cn = Pn + iQn, dn = Pn - iQu- (H) 

to the form 



'1 - 3gu(f>l - gi2'lpl)Un - 25'l20n^nPn - e(Mn+l + ^n-l) 

^1 - 5-110.^ - 5'l2^^)Wn - e{Wn+l + Wn^l) = Am„ 

^1 - 3g22i^l - gi2(t>l)Pn - ^guCpn'ipnUn " e(Pn+l + Pn^l) ■ 

[1 - g22^l - gi24>l)qn - e(gn+i + qn-l) = XPn- 



-XWr. 



-Mn 



(12) 



The matrix- vector form of the problem is then of the form: 
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£+u = — Av 
£_u = Av, 



(13) 



where 























{PnJ 





(14) 



and£± are infinite-dimensional symmetric matrices which have tridiagonal 
blocks as follows: 



(•^+)n,n 



y -2^120nV'n 1 - 3^22^^ " ^120^ 



1 - ^22^^ - C/120^ 



(15) 



±)n,n+l — y^±)n+l,n 



-e 

-e 



Equivalently, the stability problem ([T3|l is rewritten in the Hamiltonian form 

JHt7 = Xt], (16) 

where T7 is the infinite-dimensional eigenvector, which consists of 4-blocks of 
(u„,p„, w„, g„)-^, J' is the infinite-dimensional skew-symmetric matrix, which 
consists of 4-by-4 blocks of 



Jr 



n,m 



1 0^ 
-10 
1 
-1 Oy 



■'n,mi 



(17) 



and "H is the infinite-dimensional symmetric matrix, which consists of 2-by-2 
blocks of 



(£-), 



(18) 
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The representation (fT6|) follows from the Hamiltonian structure of the DNLS 
equation ([1]), where J is the symplectic operator and 'K is the linearized 
Hamiltonian. Using our above existence results, the matrix "H can be expanded 
(again for e small) into the power series 



k=l 



where 



(19) 
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(20) 



The solution of the difference equations is defined by the power series @ and 
and the pair (l)^\ipl^^ solves the inhomogeneous problem 



.,(0) 



(21) 



The symmetric matrix is defined by flTSl) . where 7/'-°^ is given by fl2Ul) and 
"H*-^^ consists of blocks: 




(22) 



where 



1,-'^+ Jn.n 



(23) 
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-2(722V^(°Vi')-2^7i20f0(^) 



-1/(1) _ nj(-L) 
'''n,n+l ''-n+l,n 



(1) 



while all other blocks of Hnln are zero. 



-10 ^ 
0-100 
0-10 
-ly 



(24) 



(25) 



Now, the non-self-adjoint eigenvalue problem f lT3|l can be transformed to the 
self-adjoint problem 

A^u = £+£_u = -A^u. (26) 



According to fl20l) . 



/r(0)^ 

)n,n 




Therefore, 



(27) 



f28) 



(29) 

where is given by (127|) and £L^'' is given by (12^ . This expansion is a 
natural generalization of the one-component results of [27]. As a result of it, 
if we identify the eigenvalues 7 of we can then trace the eigenvalues of 

the full Hamiltonian problem according to = — £7. 

We now turn to a number of special case examples, which allow us to consider 
the stability of the different two-site and three-site configurations. 



3 Case examples 



We now focus on the set S: N = 2 and = 3. Since the phases of the two 
components assume the values ^n, (n = {0, vr} for n e 5" = {1, 2, • • • , A^}, we 
use the symbolism for 0, "-" for vr and "0" for the case where no excitation 
is present on a particular site. 
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Fig. 1. The top panel shows the spatial profile of the mode ++,++ in both compo- 
nents and the corresponding spectral plane of the linear stability problem is shown 
in the middle panel for e = 0.05. The bottom subplot shows the continuation of the 
branch from e = to e = 0.15 and the real positive eigenvalue numerically obtained 
(solid line) and theoretically predicted (dashed line). The theory yields 
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When N = 2, the eigenvalue problem can be explicitly expressed in terms of 
911, 912 and (722 in the eigenvalues of Eq. f l29|) . More specifically, we obtain the 
following results. 



-+: The nonzero eigenvalues of A^^^^ are given explicitly as 
4(5'ii -5-12) (6-12 -922) 



71 



-5'12 + 5'll5'22 



72 = -4. 

Therefore, the ++,++ mode has two real unstable eigenvalues 



Ai = 
\2 = 2^fe 



\ 



4(^11 - 5^12) (5'i2 - 5-22) 



5*12 ~ 5'll5'22 



(30) 
(31) 



(32) 
(33) 



in the stability problem f|T6|) for small e > 0. Fig [T] shows the result under 
parameters gu = 922 = ^ and gi2 = 2/3, where the agreement with theory is 
very good for e < 0.05. We have generally observed this kind of agreement 
with the analytical results presented herein for an interval of e values in 
the vicinity of the anti-continuum limit. As a final observation, associated 
with this case, we should point out that although the relevant eigenvalue 
expression has non-trivial zero crossings, the constraints for the existence 
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Fig. 2. The top panel shows the mode ++,+- and the middle panel its spectral 
plane for e = 0.05. The bottom subplot shows the real and imaginary parts of the 
eigenvalues (again numerics in solid vs. theory in dashed). The theoretical prediction 



gives = ±4?e- 



of a real solution mandate that the relevant eigenvalue always stay real and 
hence the configuration be unstable 0- 
++,+-: In this case, we can obtain that 



2 



"^jgiigu - 912922) 

9 

9l2 ~ 9u922 



\/4(^ll - ^12) (^'12 - 922) {912 - 911922) + {911912 - 9l2922f 

2 ■ \ ) 

9i2 — 911922 

In this setting, an interesting possibility arises, which deviates from what 
is known about the one-component cases. In particular, the nature of the 
instability of such configurations depends on the choice of 911, 912 and 
922- When = (722 and 912 > 911, there is a quartet of complex conjugate 
eigenvalues which gives rise to an oscillatory instability, as soon as e becomes 
nonzero, shown in Fig [3]). On the other hand, when gn = (722 and gi2 < gu, 



Recall that in the Hamiltonian system under consideration here, when A is an 
eigenvalue, so is —A, A* and —A*, hence the existence of a real eigenvalue pair is 
always associated with an exponential instability. 
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Fig. 3. The top panels show the mode ++,+- but with gu = §22 = 1, 512 = 1-5 
and the middle panel its spectral plane for e = 0.05. The bottom subplot shows 
the real and imaginary part of the quartet, for which the theoretical prediction is 

a pair of real eigenvalues and a pair of purely imaginary eigenvalues are 
generated as can be seen in the case example of the continuation shown in 
Fig ED. 

• +-,+-: The mode has two purely imaginary eigenvalues 



4(fl'ii - 5'i2)(5'i2 - fl'22) 



912 ~ fi'iifl'22 



(35) 
(36) 



This is once again mandated by the existence constraints considered pre- 
viously. As a result, similarly to the one-component analog of this out-of- 
phase mode, in this two-component installment too, we find that this mode 
is linearly stable in the vicinity of the ant i- continuum limit. This mode is 
only destabilized for sufficiently large values of e ~ 0.147, when one of its 
imaginary eigenvalues collides with the continuous spectrum and leads to a 
complex eigenvalue quartet. The relevant detailed results for the continua- 
tion of this branch are shown in Fig. |H 

++,+0 and -I— ,-t-O: When [/„ is an excited site but Vn is non-excited for 
some n G S*, we have 
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Fig. 4. The top panel shows the mode +-,+- and the spectral plane for e = 0.05. 
The bottom subplot shows the imaginary part of the eigenvalues, for which the 
theoretical analysis yields = — 4e, — |e. 
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Fig. 5. The top panel shows the mode ++,+0 and the spectral plane for e = 0.05. 
The bottom subplot shows the real and imaginary parts of the eigenvalues given 
(theoretically) by A = ^i, 1.74543^6. 
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Fig. 6. The top panel shows the mode +-,+0 and the spectral plane for e = 0.05. The 
bottom subplot shows the two pairs of purely imaginary eigenvalues theoretically 
predicted as A = 1.74543iVe, They collide from e = 0.025 to e = 0.041. 
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so that 7^ 0, which induces an nonzero eigenvalue A = — or 
|— — respectively, in the eigenvalue problem f[T^ of order e*^. This 

922 

already provides the dominant behavior of the corresponding eigenvalues 
as can be observed in the eigenvalue plots of Figs. and O In addition 
to this purely imaginary eigenvalue, there is a small eigenvalue bifurcating 
from the origin in both of these configurations, which is well approximated 
(for gu = 922 = ^ and gi2 = 2/3 in the figures), by A^ = ±1.74543e. 
However, this higher order eigenvalue plays a critical role in the stability of 
the configuration, as in the case of the ++,+0 configuration this eigenvalue 
exists as real (i.e., with the + sign in the above expression), while in the 
latter case of +-,+0, it is imaginary. Thus, the former is unstable and the 
latter is stable, at least for small values of e; for larger e, as can be seen in 
Fig. [6l a collision arises between the ©(e") near-constant eigenvalue and the 
growing O(v^) eigenvalue (between e = 0.025 and e = 0.041) which gives 
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Fig. 7. The mode +++,+++ (top panel) with four pairs of real eigenvalues (see the 
spectral plane of the middle panel) is shown in the figure, captured when e = 0.05. 
The analytical prediction for these pairs = 6e, 2e, |e, |e is shown by dashed 
line in the bottom panel and compared to the numerical results (solid lines). 

rise to a complex eigenvalue quartet and hence an oscillatory instability. 

When = 3, analytical expressions for the eigenvalues A are too cumbersome 
to obtain/write explicitly, however the relevant eigenvalues can be computed 
(as a function of e) by solving a low dimensional eigenvalue problem for each 
particular set of Qij. Fig [71fT3] summarize the results for the three modes under 
the parameters gu = §22 = ^ and gi2 = 2/3. 

The case of the +++,+++ configuration is shown in Fig. [TJ where we can 
again verify the relevant expectations: for each excited site, we expect one 
pair of emerging eigenvalues, hence 6 pairs should be associated with the 3 
excited sites in each component [notice that in the case of 2 excited sites, this 
amounted to 4 such pairs]. Out of these, 2 are always set at the origin due to 
the phase invariance of each of the involved components, hence 4 pairs remain 
as non- vanishing [2 in the 2-excited site case]. The in-phase structure of the 
configuration renders all of them real in this case in good agreement with the 
predicted theoretical eigenvalues (of A^ = 6e, 2e, |e, |e given by the dashed 
lines in the graph). 

The other two configurations which are identical across components are the 
++-,++- of Fig. [8] and the +-+,+-+ of Fig. |9l In direct analogy to their one- 
component counterparts examined in [27], we find that the former of these two 
configurations bears 2 real pairs of eigenvalues and 2 imaginary ones, while 
the latter is indeed the only 3-site configuration that enjoys linear stability 
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Fig. 8. Similarly to Fig [7] but for the mode ++-,++- with two pairs of real and two 
pairs of purely imaginary eigenvalues given by = ±2y/3e, ±^e. The top three 
panels are captured when e = 0.05. 




Fig. 9. Similarly to Fig [7] but for the mode +-+,+-+ with four pairs of purely 
imaginary eigenvalues theoretically given by = — Ge, — 2e, — |e, — |e- The top 
three panels are captured when e = 0.05. 
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Fig. 10. The left column shows the mode +++,++-, which terminates at e = 0.072 
by colliding with the mode +0+,++- shown in the right column. The top two rows 
show the profiles of each mode when e = 0.05. The bottom two rows show the 
eigenvalues of each mode. The mode +++,++- has three pairs of real and one pair 

of purely imaginary eigenvalues given by = I (6± \/2l) e, ±2Y^e, while the 
mode +0+,++- has two pairs of real and two pair of purely imaginary eigenvalues 
givenbyA = |i, 1.61812iY^, 1.61812y^, 1.82599e. The first two pairs of eigenvalues 
collide and form a complex eigenvalue quartet from e = 0.025 to 0.048. 



in the vicinity of the anti-continuum limit, due to the out-of-phase nature of 
the neighboring sites therein. In the case of the ++-,++- configuration, the 
eigenvalues are analytically found to be: = ±2-\/3e, ±^^e, while in the 
case of the +-+,+-+, the A^'s assume the opposite value than they do for the 
+++,+++ case (yielding 4 imaginary pairs). 

Finally, we consider configurations which are distinct in their 3-site form be- 
tween the two components. These are shown in Figs. ITOlfTSi It should be noted 
that some of these configurations cease to exist within the interval of e (up 
to 0.15) probed herein. In particular, for instance, the mode +++,++- of 
Fig. [To] collides and disappears with the branch +0+,++-, as can be seen 
in the figure; nevertheless in its interval of existence, we can obtain good 
agreement with the theoretical prediction for the corresponding eigenvalues 
A^ = I ± \/2T) e, ±2y^e. In turn, for the mode +++,+-+, we find that 
it collides with +0+,+-+ around e = 0.06; the details of the corresponding 
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Fig. 11. Similar as Fig. [10] but here the left column shows the mode +++,+-+, 
which has eigenvalues A'^ = ±-^e, i"^^ the right column contains the mode 

+0+,H — h, which has eigenvalues A = ^i, 1. 24467^6, 2.14026i\/e, 1.75819e (ac- 
cording to the theoretical analysis). The first two pairs of eigenvalues collide and 
form an eigenvalue quartet from e = 0.014 to 0.028. The top four panels are captured 
when e = 0.05. 

eigenvalues are shown in Fig. [11] Lastly, the mode ++-,+-+ collides with 0+- 
,+-+, as shown in Fig. [121 and the mode ++-,H — collides and disappears with 
0+-,+-0, as detailed in Fig. [T3] 



++-,H — becomes 0+-,+-0 after e = 0.1, etc. 



4 Conclusions &; future challenges 



In the present work, we considered two-component settings and provided a 
general framework for the existence and stability of localized mode solutions 
to nonlinear dynamical lattices of the DNLS type. We illustrated that and 
TT relative phase configurations arise in the context, in analogy to the case of 
the one-component DNLS model. Although the stability presented some inter- 
esting deviations from the corresponding one-component counterpart (e.g. in 



17 




0.5 



0.5r 



0.04 0.08 0.12 

e 



0.04 0.08 0.12 

£ 




"0 0.04 0.08 0.12 

e 



0.04 0.08 0.12 

e 



Fig. 12. Similar as Fig. [10] but here the left column has the mode ++-,+-+, which 
has eigenvalues = — | (6 it \/2T) e, while the right column contains 

the mode 0+-,+-+, which has four pairs of purely imaginary eigenvalues given by 
A = \i, 0.834064Ve, 1.39576iVe, 2.27766iv/e. The eigenvalues of order e° colhde 
with the others from e = 0.016 to 0.023, from 0.039 to 0.052, and from 0.063 to 
0.102, respectively, giving rise to the corresponding oscillatory instabilities. The top 
four panels are captured when e = 0.05. 



the case of the ++,+- configuration), most of the cases presented significant 
analogies in the nature of the stability conclusions (with suitable modifications 
of the eigenvalue counts) to the one-component setting. In particular, the ex- 
istence of in-phase near-neighbor excitations in the same component appeared 
to generically lead to the existence of real eigenvalue pairs in the linearization. 
On the other hand, the only configurations that seemed to generically be sta- 
ble were those where all adjacent sites were out of phase with respect to each 
other. 



It is natural to envision numerous generalizations of the present considera- 
tions. While already some special case examples of discrete vortex configura- 
tions have been examined in [T7], a general theory of both discrete solitons 
and of vortices is still absent in the two-component setting. Far less work has 
been done in the three-dimensional case, where identifying suitable stable con- 
figurations with nontrivial phase profiles is a fairly challenging task even in 
the one-component case [19]. While these considerations are still lacking for 
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Fig. 13. Similar as Fig. [10] but the left column shows the mode ++-,H — , which has 
four pairs of eigenvalues = ±2-y/|e, :ii2^J^e, while the right column contains 
the mode 0+-,+-0, which has four pairs of purely imaginary eigenvalues given by 
A = -^i, 1.44578i-v/e, 2.00338«-y/e. The eigenvalues of order collide with the 
others from e = 0.018 to 0.039, from 0.036 to 0.051, respectively, creating once again 
an oscillatory instability. The top four panels are captured when e = 0.05. 

the more standard DNLS model, this is certainly all the more so the case for 
the setting with four- wave mixing terms being present, as appears to be ex- 
perimentally relevant per the work of [20j . Some of these avenues are presently 
being considered and will be reported in future publications. 
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